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It is shown that the fundamental standard hadron cascade diffusion equation in the Mellin trans- 
form space is not rigorously correct because of the inconsistent double energy integral evaluation 
which generates the function < n s > with its associated parametizations. To ensure an exact basic 
working equation, the Faltung integral representation is introduced which has the elasticity distri- 
bution function u(rj) as the only fundamental input function and < r/ 3 > is just the Mellin transform 
of u(rj). This Faltung representation eliminates standard phenomenological parameters which serve 
only to mislead the physics of cascade. The exact flux transform equation is solved by the method 
of characteristics, and the hadron flux in real space is obtained by the inverse transform in terms 
of the simple and essential residues. Since the essential residues are given by the singularities in 
the elasticity distribution and particle production transforms that appear in the exponentials, these 
functions should not be parametized arbitrarily to avoid introducing non-physical essential residues. 
This Faltung formulation is applied to the charged hadron integrated energy spectra of halo events 
detected at Mt. Chacaltaya by the Brazil- Japan Collaboration, where the incoming flux at the 
atmospheric boundary is a single energetic nucleon. 
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1. INTRODUCTION 

The Brazil- Japan Collaboration detected many cosmic-ray events in the energy range (10 13 — 10 17 ) eV with emulsion 
chamber exposed at Mt. Chacaltaya. About twenty events were observed in the visible energy region E > 1000 TeV. 
Approximately, half of them p], |2, H, E| were associated with a uniformly darkened wide area on X-ray films. The 
central part of this area was called " halo" , and so these events were called " halo events" . Similar events were also 
observed in experiments at Pamir [5[,Fuji and Kanbala Q- Recently, a new experiment Q, using a hadron 
calorimeter associated with emulsion chambers at Mt. Chacaltaya, reported this kind of events as well. Thus, the 
appearance of a strong concentration of energy and particles as a halo seemed to be a common feature in this energy 
region. After the 1980s, these super-families were also compared with simulations using different primary compositions 
and models for high-energy nuclear interactions 0, EH EH- Nevertheless, these simulations could not describe fully 
all the events with the same inputs on primary composition and nuclear collision models. Due to their large hadronic 
number and energy, some authors [l2| suggested that these events could be explained as centauro-like ones. 

Recently, we have developed an analytic method [l3| that allows us to calculate the hadronic and electromagnetic 
components of cosmic rays in the Earth's atmosphere, and generalizes earlier works [T3, EH EH E3] by including the 
energy dependence of the mean- free path effect. An important issue in the high energy region concerns the behavior of 
the inelasticity, which is defined as the fraction of energy given up by the leading hadron in a collision induced by an 
incident hadron on a target nuclcon or nucleus. This parameter has been exhaustively studied in several papers, but 
until now it continues to be an open question. Several authors have suggested that the average inelasticity coefficient 
is an increasing function of the energy , whereas others have proposed that it is a decreasing one 

High energy cosmic rays, which reflect the nuclear interaction in the energy region covering 1 to 100 TeV, are well 
fitted with a constant mean inelasticity 1/2 [23l. 1 241] . At higher energies, where these super-families belong, a constant 
value for the inelasticity is no longer valid in order to explain experimental data. Up to now, all papers treating 
analytically this subject have used simplified inputs on interaction mean-free paths and inelasticities. 

Here, we believe much of the difficulties in developing an analytic approach of cascade lies on the fact that the 
widely used standard diffusion equation in Merlin transform space is not a mathematically rigorous equation, due to 
the approximations used in the evaluation of the double energy integral that brings in the entity < if >. Because 
of this fundamental inconsistency at the very beginning, it is reasonable to understand that each level of the cascade 
would compound the array of free parameters. In section [2] we make use of the Faltung theorem to formulate exactly 
the diffusion equation in Mellin transform space. In this exact Faltung representation, < 77 s > is just the transform 
of the elasticity distribution function u(rf). There is no need to model and parametize < r/ s > itself. Consequently, 
elasticity distribution is the primary and only input function. As a result, we could avoid over parametizing and 
contaminating a model based on an approximated non-rigorous transform equation. Method of characteristics (l3j 
is used to solve the exact transform equation considering one single incoming nucleon as the boundary condition in 
section [3J For the nucleon case, the flux in real space is evaluated by the s = — 1 essential residue in section SJ Our 
solutions are presented in the usual modified Bessel functions of order 1. For the pion case, the Faltung representation 
of the transformed equation is treated in section [5j and the s = sq simple and s = and s = — 1 essential residues are 
evaluated in section^ In sectional we compare these fluxes with the halo event data of P06, Ursa Maior, Andromeda 
and Mini- Andromeda III @, d, Sp discussions and conclusions. 
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2. FALTUNG FORMULATION OF NUCLEON DIFFUSION 



From considerations of different fundamental physical processes, the number density flux energy distribution func- 
tion of nucleons N(E,t) at energy E and atmospheric depth t is described by 

■ r U (E/ E ')^4^ dE ' (i) 



HE) J E HE') E> 

where X(E) is the energy dependent mean-free path, r\ = E/E' < 1 is the elasticity, u(E/E') is the elasticity 
distribution. Modelling the mean- free path by a power index (3 [25(, 

\{E) = \ N {^)-P, (2) 
the second term on the right side of Eq. (fT]) can be converted from a E' integral to an ry integral that reads 



dN(E,t) 

Ft 



\ N B X N B J Q r\ r\ 



(3) 



where B is the normalization energy of the mean-free path. We note that the energy dependence in X(E) makes the 
mean free path decrease as energy increases, (E/B) > 1. The power index (3 > depends on the reference energy B. 
We observe that the nucleon cascade equation, Eq. fl}, nas twocompeting terms on the right side. The first term is 
the diffusion term that drains the flux N(E, t) dE at E to lower energies E'. The second term is the attenuation term 
that fills the flux at E by higher energies E'. Since the mean-free path scaled by Eq. @ vanishes as E/B goes to 
infinity with 0, the first term would dominate the equation and the spatial gradient of the flux would be very negative 
at high energies. 

Instead of solving Eq. ([3]) in real space [HI, we use the integral transform approach. Before doing the 77 integral, 
we take the Mellin transform defined by 

N(s,t) = j t)d(-) (4) 

where the energy E is normalized to some reference energy A, so that the transform does not carry dimension of 
energy to power s. Now, Eq. ([3]) in the transform space reads 

at Xn B 

Xn B J Q 77 77 77 

To reach this equation, we have done the 77 integral locally to get < rj s > while ignoring the 77 dependent terms in the 
E integral. With K as the normalization energy the following model of parametization is often used (2lj | 

< v s >= / u{ri)rfdn = n ■ . , (—) KS 



1 1 A E' A F' 

:( r( r = a(s)( r( )KS (?) 



(s - s ) S y K' A v ,X K' K A 

where s = —1/8. For a uniform elasticity distribution, we have k = 0, 5 = 1, and s = — 1- I n particular, taking 
s = 1 gives the average elasticity 
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Since < r/ s > is parametized to an E' dependence, the E integral is now converted back to an E' integral with the 
help of r\ to complete the second term. The equation of the flux transform then becomes 

Ajy K 

We notice that, should the energy E in the Mellin transform be not normalized to some reference energy A, then 
N(s, t), N(s+(3, t) and N(s+(3+ks, t) would have different dimensions in energy. Here, in Eq. ((BJ), they have the same 
dimension of N(E, 0). We remark that the second term on the right side of Eq. (J6j is obtained by isolating the < rj > 
integral to get < 77 s > and leaving other rj dependent terms to be converted to the E 1 integral. Such mathematically 
unwaranted procedure makes Eq. ([6]) not an exact representation of Eq. ([1]). Consequently, any parametization on 
the elasticity in the approximated flux transform equation could only contaminate the physics of the nucleon flux in 
real space. 

To avoid this serious mathematical inconsistency, we reconsider Eq. |T]) and make use of the Faltung integral. By 
defining the elasticity distribution u(E/E') = for E/E' > f , the lower limit of the integral can be extended to zero, 
and Eq. ([T]) reads 

dN{E/A,t) N(E/A, t) [°° JV(g /A,t) 1 , 

f°° 1 

= -w(E,t)+ / u(E/E')w{E',t)—dE' 
Jo E 

= -w(E,t)+u®w(E,t) (10) 

The last equality is reached by using the Faltung integral of the Mellin transform. Taking the Mellin transform leads 
to 

9N(s,t) _, . , „, \ 

= — w(s, t) + u{s)w{s, t) 



Of 



-[1 - u(s)]w(s, t) = - J_(£y [1 - u(s)}N(s + 13, t) (11) 
an £> 



poo 771 rp rp 

Hi*) = J o (-y u (-)d(-)=< v s > (12) 

By changing the normalization energy A in Eq. (fl2]l to E', and recalling u(E/E') = for E/E' > 1, we see that 
< 77 s > is just the Mellin transform of the elasticity distribution function. It is completely determined once elasticity 
distribution is specified. 

In this Faltung representation, elasticity distribution 11(77) i s t ne fundamental parameter. There is no need to 
specify < r/ s > together with those unnecessary parameters associated to it. The necessity of parametizing < 77 s > 
is actually derived from the approximated nature of Eq. This could over parametize the model based on an 
approximated equation, which could generate inconsistencies among other parameters. An exact representation of 
Eq. (TTJ) by Eq. (jlip through Faltung theorem sets a firm base upon which additional features of one-dimensional 
cascade could be analyzed with confidence. This Faltung representation apparently requires u(E/E') as a function 
of E/E' only. However, this is not true. With rj and E' dependences in a separable form, we could transfer the E' 
dependence of u(r),E') to w(E',t). By defining a new w(E',t), the Faltung formulation prevails. 
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3. METHOD OF CHARACTERISTICS 

With the Faltung formulation, we proceed to solve Eq. (fTTj) . Some researchers solve equation of this kind formally 
by operators [13] ■ Since (3 is much less than s, we choose to make a Taylor expansion of N(s + (3, t) to get a first order 
differential equation 

X N (^^ + [l-u(s)]0 d -^ = (13) 

This partial differential equation is equivalent to the following set of ordinary differential equations which describes 
the trajectory of the coordinate point (s, t, N) in the functional space parameterized to £ (28l. l29l [30| 

(±)P*1 = ds = dN 

V \ N [l-u(s)}(3 [l-u( s )]N 1 ' 

This method of characteristics to solve first order partial differential equations was used in superradiant free electron 
lasers [29|, [3(J. Solving for the equality between dt and ds, 

( A dt_ ds 

K B> X N [l-u(s)]P [ ' 

we get a trajectory between the variables t and s through the parameter £, t = t(s, (3, k), which is the characteristics 
of the partial differential equation, Eq. (flU)) . 

To get the transform of the flux, we could solve the equality of dN with d£, or with ds, or with dt. Since the 
boundary condition of N is given in terms of s at t = 0, we choose to solve with dt 

A f3 dt _ dN 



B' \ N [l-u(s)]N 

N{s, t) = N{s, 0) e -{A/Bf[i-u{ s ))t/x N (16) 

For the boundary condition, we use a single incident nucleoli with energy Eq as follows plus its associated Mellin 
transform 

*<f'°> (17) 

N(s,0) =(^r s . (18) 
Eq 

The inverse transform of Eq. (|16[) . therefore, gives the nucleon flux in real space 

JV(^ >t ) = (^°)-lJ_ f(]L)-{s+l) e -{A/Bni-u{s)]t/X Nds _ (19) 
A A 2ni J Eo 

For the elasticity distribution, which is our fundamental parameter, let us consider a linear profile 

u(rj) = 1. (20) 

By taking the Mellin transform, we get 

u(s) = J 1 I . (21) 

We note that the parameter j3 which measures the mean free path effects appears in the (A/B) 13 factor and in the 
characteristics of the partial differential equation. 



6 



4. CALCULUS OF RESIDUES 

To get the nucleon flux in real space, we need to evaluate the residues in Eq. (|19p . An inspection of the equation 
tells that there are no simple poles. However, in the exponent, there is a pole at s = — 1 in u(s) 

u(s) = -L- (22) 
s + 1 

Since it appears in the exponent, this is an essential singularity. The only contribution to the nucleon flux, therefore, 
comes from the s = so — — 1 essential singularity. To consider this essential residue, we expand the exponential function 
in power series to obtain 

U — W P -(A/Bft/X N 



2mRes(s ) = (-i)-i {(_)-* e 



(23) 



We define the function acting on the (s — sq) powers inside the integral by G(s) = (Eo/E) s = g s . Since G(s) is 
analytic in the neighborhood of s — so, we expand it in a Laurent series about s — sq so that Res(so) becomes 

2mRes(a ) = [E)^ {e -{A/B)H/x N 

y y ±_ J_ * ]n G M ( ) f — i — ds} _ (24) 

By taking (n — m) = +1 terms, we pick up the contributions to the essential residue so that 

n— 1 v 7 

Here Z 2 = 4(A/B)P(t/\ N )(\ng). The infinite series can be rewritten in terms of the modified Bcsscl function of order 
one. Thus, the result is 

w( f''' = (|)- I (J)-"«- < " /B) "" A " 

We recall that the series that gives rise to the modified Bessel function of order one is a semi-divergent series. It 
diverges up to some nth term due to the factors like (t/Ajv) and (hig). Afterwards, it begins to converge due to 
the factorials of n. For the single incident nucleon case, the fraction (Eq/E) > 1 so that g is larger than unity and 
(In g) is positive. Consequently, Z 2 (n) is positive which leads to the modified Bessel function solution. The essential 
residue here plays a very important role, because it represents the flux at a given atmospheric depth i/Ajv. We remark 
that this case of single incident nucleon had been solved in real space in terms of probability distributions under the 
assumption of constant mean- free path and uniform elasticity [lj, [l5| . It was also solved with Mellin transform by 
residues under the same assumption [l7j |. In both approaches, the results are in terms of the modified Bessel functions 
h{Z). To obtain the integrated flux, we use 



E i-Eo/A E E 
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5. PION DIFFUSION 

For the pions, the number density flux energy distribution function U(E,t) at energy E and atmospheric depth t 
is given by 



dU{E,t) U(E, t) 

at ~ ~ \'(e) 



+(i - t >r a(E/E ' ) iw^ E ' 







vAE/E) l^YW dE + J Vn{E/E) ^(EYW dE (27) 



Here b = 1/3 is the charge exchange probability of the incident pion, v(E/E')/E' — <j)(E' , E), and <f>{E' , E)dE is the 
number of produced particles at energy E in the interval dE due to the incident particle at energy E' . Furthermore, we 
have again defined u{E/E') = and v(E/E') = for rj = E/E' > 1. We also have assummed that v(E', E) = v(E/E'). 
This allows us to use the same Faltung formulation to write 

d 11 ^' - = -w n (E, t) + (l- b)u ® w n (E, t)+ V7I (g> w n (E, t)+v n ® w(E, t) (28) 
The Mellin transform of the pion equation reads 

d n |^' - = -w w (s,t) + (1 - b)u(s)w v (s,t) + v n (s)w n (s,t) + v n (s)w(s,t) 
at 

= -[l-(l-b)u(s)-v 7T (s)]w(s,t)+v n (s)w(s,t) (29) 
Again, we model the mean-free path of pions by the same power index [3 

A'(£) = A 7r (|)-' 3 (30) 
As for the pion and nucleon production by incident pion, there is the Feymann scaling which reads 

4>(E',E) = D[l-(^)] d ± (31) 

where D = (d + l)/3, d = 4. This gives <j>(E', E)E' = v(E/E') which allows the application of the Faltung theorem. 
Using the method of characteristics with the boundary condition n(s,0) = 0, the pion flux transform is 

P{s)/K - M(s + P)/\n 

X {1 - e +(A/B)^(s+f3)t/\ Ne ~(A/BfP(s)t/\^ ^ 

ft(s) = [1 - u(s)} (33) 
P(s) - [1 - (1 - b)u(s) - Q(s)] (34) 

Q(s) = v(s) (35) 

This pion transform is expressed in terms of the nucleon transform which acts as a source term. Here, we note that 
(3 appears explicitly in /i(s + (3) and N(s + (3,t). The mean free path effects are explicit in the pion flux. Doing the 
inverse transform gives the pion flux in real space 
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m E t) = -^-A-i [( E r (s+1) q{s)/Xn 

K A' > 2ni y A> J y Eo> -P(s)/K + n(a + 0)/X N 

x {_ e -(A/B)^{s+l3)t/\ N + e -(A/BfP(a)t/X„y ds 

- r?(T ) _1 A#)~ (s+1) ^) 

Zm A J Eq 

x r_ e -(A/B)f>rts+i3)t/\ N +e -{A/B)f>p(s)t/K\ fa (36) 

As has been discussed in an earlier publication [l3| , the pion flux is given by the residues of the simple singularities in 
Z(s, 0) and the residues of the essential singularities in the exponents /x(s + (3) and P(s). These essential singularities 
are expressed in an infinite series with terms involving higher and higher derivatives. For the nucleon case, there 
happens to be a closed analytic form for the essential residues as in Eq. ([26]) . In most cases, we have to evaluate them 
numerically 

Since the essential residues depend on the singularities of <fc(E', E) in the exponent, the functional form of </>(£", E) 
acquires a real physical importance, beyond simple data fitting. For example, should we choose Eq. (|3ip . we would 
have v{rj) = <j>(E', E)E' = (5/3) (1 - ry) 4 /?? and 

Q(s) = v(a) - D ^(-lrC*^— (37) 

There would be five essential singularities s=0, -1,-2, -3,-4, and each singularity would generate a pion flux represented 
by its residue. Some of the residues are extraordinarily large and others are extraordinarily small. In our opinion, 
these residues are due to singularities that are introduced by data fitting, whose fluxes in real space are unphysical. 
We could also choose 

v(n) = 2.08-(l -ii)e- b11 
V 

This would generate a Gamma function dependence in v(s) which would give an infinite sequence of singularities 
on the left side of the complex s plane. Consequently, under the method of residues, arbitrary fitting of v(s) could 
contaminate the flux. 

In order to bring out the physics, we consider the energy conservation between the average leader pion and the 
average secondary pion energies 



<r]> + <K> = J [rfu(rj) + rjv{rfj\dri = 1 (38) 
To obtain the distributions u{vj) and v(rj), we note that the integrant should satisfy 

r}[u(rj) +v(rj)] = l,2?7,3ry 2 , ... 
Taking the first choice on the right side, a consistent pair of distributions is 

u(ri) = (1 - arf) (39) 



v(rf) = -(1 -77 + arj 2 ) (40) 



with elasticity and inelasticity given by respectively 



// = / r]u{ri)drj = \~\ (41) 
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<K>= J Tjv^drj = \ + l (42) 
Taking the Mellin transform of u(rj) and v{rf) gives 

u(s) = [t-^— - t-^] (43) 



This energy conservation consideration identifies two singularities of v(s) at s=0,-l with a = 0. Again, the assumption 
that v(E\E) — v(E/E') appears to be restrictive. The energy dependence of v can be introduced in the Faltung 
representation in the same way as u which has been discussed earlier. 
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6. PION RESIDUES 

To evaluate the pion flux of Eq. (j36|) , we have to consider the simple residues of Z(s, (3) and the essential rresidues 
of fi(s + 13) and P(s) in the exponents. The nature of these residues have been discussed before (l3j . For the present 
case with u(s) and v(s) given by Eq. and Eq. (|4"4")) respectively with a = 0, using a — X^/X v = 5/7, and writing 
Z(s, (3) = Z(s) by neglecting the mean free path effects in this algebric function, the simple poles of Z(s) with 6=1/3 
are given by a pair of complex conjugate roots 

(so, s* ) = (s r ± isA = ^(5 ± z^) 1 / 2 ) (45) 

As for the exponential function, the mean free path effects are represented by jj,(s + (3) in the first exponent of Eq. (I3GD . 
Here, even j3 is retained to evaluate the exponential, the effects are still small on a Log-Log plot. Evaluating fj,(s) and 
P(s) at these roots give respectively 

(/V ± ifn) = 1 - | — -|2((s r + 1) =F«s*) (46) 

I so + ir 

(P r ± iP*) = 1 — |— j2 ( s r T i*<) + _l_(( ar + 1) T i Si ) (47) 

kor |so + i|- 

With these notations, the simple residues are 

So ,_i rJ E , (sr+1) 1 1 



2^(i?e S ( So )+i?e S (^)) = (-^)- 1 ( — ) 



A y E ' {I -a) Si 

x {e -MA/Bf mN) sin[fl A f _L + s . H E }] 

- e -Pr(A/Bf(t/\«) R - n{p jA^ t , E 



sin[P i (-f-+ Si ln(— )]} (48) 



These simple residues have a (E/Eo)~( Sr+1 ' ) power dependence on energy, plus a ln(E/E ) energy dependence in the 
argument of the sine function. Because of these dependences, this part of the differential pion flux on a Log-Log plot 
is a straight line with — (s r + 1) slope plus a slight modulation due to the sine function. 
As for the essential residues, let us rewrite Eq. ([56")) as 

n(4,t) = -^(-2)-i e -w^('A«) / ( _^)-(s+i) z(s ) e +«(s+w*/A«) ds 

A 2ni A J Eq 

+ J_iE^\-i e -(A/B)f , (t/\,) I r^y{s+X) z r s y+[{i-b)u(s)+v( s )]{t/K) ds 
2tti A J E 

Substituting the expressions of u(s) and v(s), we get 

n(^, t ) = _J^)-l e -(W(tA*) [ { ]l r (s+l) Zis)e HA/Bf(l/(s + +m t/X N ) ds 

A 2tti A J Eq 

+ _L ( ^ ) -i e -(^)^ ( ^)| ( |. ) -( S+ i) z(s) 

xe +{A/Bf(+l/s)(t/K) e +{A/B)f , {-b/( s +mt/K) ds 

= -^-(^)- 1 e - (A/B)3(tA " ) J(§r^Z(s)W 1N (s,x)ds 

+ ^(^r 1 ^ (A/B) " ( * /A - ) /(^)-( s+1 ^( S )H/ .(s, a; )^( S ,a ; )d S (49) 

Z7TZ A J Eq 

The exponent associated to in the first term has a pole at (s + (3 + 1) = (s + 1) =0. Following the same 
mathematical procedures of nucleon flux of Eq. (|25p . we have the partial flux 
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2mResN{-\) = -(Eyi e -(A/Bf(t/x N ) 

00 1 1 A t 

x{E^^Kb)^]™ + 'G«(»,,)> (50) 

G 0N (s,x) = (^yZ(s) (51) 
The two exponents associated to in the second term have poles at (s + 1) = and s = 0, and the partial fluxes are 

2niResn(-l) = +(1)-^-^'^^^ 



x 



<5Z h^Kb^^'oSW)} «*> 

rra=0 



2mResn(0) = +(|)-i e -(W(*/A.) 

x{Ei(^4» s i]" + ' G - <«-»»» 

m=0 v ; T 

G( 0> i) w («, a:) = (fr) s -Z(s)WW( s . 35) (54) 

Since the essential residues are derived from u(s) and v(s), the functions u and v have real physical significance. They 
directly determine the flux in real space. For this reason, u and v must not be arbitrarily fitted to data, much less in 
raising the power to avoid generating artificial residues. The integrated pion flux is given by 
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7. RESULTS AND CONCLUSIONS 

Figure 1 shows the integral hadron fluxes of four halo events detected at Mt. Chacaltaya by the Brazil- Japan 
Collaboration. In order to make numerical calculations about the integral hadron flux with our model and to compare 
with the detected events, we take the mean-free path normalization energy B = ITeV, \jy = 80 g/cm 2 and ft = 0.06 
which are obtained from accelerator and EAS data in the region 1 TeV < Ei a b < 1000 TeV (2||. For the pion mean- 
free path, we assume that \n/^tt = 5/7, and that it has the same energy dependence like the nucleon case, Eq. @ 
[3ll. [3^] . In the present calculations, we have only one free parameter, A(TeV), which is the normalization energy in 
the Mellin transform, Eq. and Eq. ([5]). This parameter is used to repreesnt the normalized energy E/A. In order 
to be compatible to the horizontal axis of Fig.l, we choose A = Eq = 1000 TeV so that E/A < 1. With all these 
parameters and considering a uniform elasticity distribution and constant mean free path with ft = 0, Fig. 2 shows the 
integrated nucleon flux labelled by N from Eq. (|26|) . the integrated pion flux from nucleon labelled by PN due to the 
s = — 1 essential residue from Eq. (|5D|) . the integrated pion flux from pion labelled by PPe due to the s = essential 
residue from Eq. (|53p . the integrated pion flux from pion labelled by PPs due to the s = sq, Sq simple residues from 
Eq. We note that the integrated pion flux from pion due to the s = — 1 essential residue from Eq. (f52"| is not 

plotted because it is much smaller than the others. Also the integrated pion flux from pion due to the s = Sq, Sq simple 
residues in Fig. 2 shows three bounces. This is not true. Actually, the second bounce corresponds to negative values, 
and PPs is an oscillating function. This is clear from the energy dependent argument of the sine function in Eq. (|48[) . 
Since we are doing a Log plot, the absolute value is taken which rectifies the negative part and gives the pulsating 
appearance. The overwhelming contribution of the integrated hadron flux comes from the PPe s = essential residue. 
The total integrated flux with ft = is shown in Fig. 3 which presents a curvature compatible to Fig.l. We note that 
Fig. 3 is the integrated flux in terms of the energy distribution function which could be multiplied by an arbitrary 
constant. This amounts to displace the plot along the vertical axis to the scale of Fig.l for comparisons. To show 
the effects of energy dependent mean free paths, the total integrated flux with ft = 0.0G is also shown in Fig. 3 for 
comparisons. 

To conclude, we have pointed out in Sec. 2 that the generally accepted diffusion equation in the Mellin transform 
space in terms of the parameter < rj s >, Eq. ^ and Eq. ijTj). is not a mathematically rigorious representation because 
of the approximations involved in evaluating the double integral in energy. Since two wrongs do not make a right, the 
subsequent parametization of < rj s > does not correct the mistake, and the model gets distorted from the start. Since 
this diffusion equation is the fundamental equation in evaluating fluxes of different generations where each generation 
relies on the flux of the prior flux of the parent generation, any misrepresentation on this basic equation would cascade 
accumulatively rendering the final results entirely meaningless. To overcome this primary problem, we have used the 
Faltung theorem to formulate an exact self-consistent diffusion equation for the flux transform, as given by Eq. (fTTj) 
and Eq. (|12p in Sec. 2. This Faltung formulation requires only the elasticity distribution, not < 77 s >, as the primary 
input function which avoids excessive, and often conflicting, parametization. Energy dependence on the elasticity 
distribution could be included in this Faltung formulation. 

The flux transform is solved by the method of characteristics, and the flux in real space is evaluated by residues, 
simple and essential. Since the essential residues come from singularities in the exponent which contains u(s) and v(s), 
this makes the choice of u(rj) and v{rf) particularly important. Arbitrary profile fitting could bring in more residues 
that might not be physical, such as Eq. (f3"T|) . Care should be exercised not to raise the power of 77 in (j>(E',E) for 
fitting purposes. This would increase the singularties in v(s) generating artificial fluxes. Although we have considered 
r\ that is independent of the incident energy E', this Faltung formulation does allow generalizations to rj(E') cases. 
This self-consistent Faltung formulation provides a firm and reliable starting point for one-dimensional cascades. 

Using a method recently developed by us [l3j], we have calculated the nucleon flux through the essential residue 
at different depths in a wide energy range initiated by one single nucleon. We have generalized earlier results to 
include the energy dependence in the collision mean-free path. Our solution is presented in the usual modified Bessel 
functions of order 1 with an energy dependent mean free path argument. For the pion flux, it is given by the simple 
and essential residues as well. A comparison of the integrated hadron flux with the hadronic spectra measured at 
Mt. Chacaltaya for four halo events shows good agreements even with ft — 0. Our hadron flux presents a curvature 
on the Log-Log plot just like the observations, althought we have used a simple uniform elasticity distribution. The 
fundamental reasons of this good agreement stem from three aspects. The first is the exact representation of the 
diffusion equation through Faltung formulation with elasticity distribution as the basic input function. This exact 
formulation eliminates unnecessary free parameters that only serve to contaminate and distort models such as < rj s > 
which is forced upon by a bad derivation. The second is the use of essential residues to evaluate the flux which has 
been overlooked by earlier investigators. The third is the use of elasticity distribution u(rf) and consequently the 
self-consistent v{rj) to evaluate the flux. Because of these three aspects, our model with no free parameters gives 
results compatible to observations using only a simple uniform elasticity distribution. 
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FIG. 1: The experimental data of four halo events at Mt. Chacaltaya. 
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FIG. 2: The (3 — integrated flux of the nucleon (N), the pion from nucleon of the essential residue s — (PN), the pion from 
pion of the simple residue (PPs), the pion from pion of the essential residue s = (PPe). 




FIG. 3: The (3 = and f3 = 0.06 total integrated pion flux. 



